Fluctuating Nematodynamics using the Stochastic Method of Lines 
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We construct Langevin equations describing the fluctuations of the tensor order parameter Qap 
in nematic liquid crystals by adding noise terms to time-dependent variational equations that follow 
from the Ginzburg-Landau-de Gennes free energy. The noise is required to preserve the symmetry 
and tracelessness of the tensor order parameter and must satisfy a fluctuation-dissipation relation 
at thermal equilibrium. We construct a noise with these properties in a basis of symmetric traceless 
matrices and show that the Langevin equations can be solved numerically in this basis using a 
stochastic version of the method of lines. The numerical method is validated by comparing equi- 
librium probability distributions, structure factors and dynamic correlations obtained from these 
numerical solutions with analytic predictions. We demonstrate excellent agreement between numer- 
ics and theory. This methodology can be applied to the study of phenomena where fluctuations in 
both the magnitude and direction of nematic order are important, as for instance in the nematic 
swarms which produce enhanced opalescence near the isotropic-nematic transition or the problem 
of nucleation of the nematic from the isotropic phase. 

PACS numbers; 05.10.Gg, 02.50.Ey, 02.60.Cb, 64.70.mf 
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I. INTRODUCTION 

Fluctuation phenomena in nematic liquid crystals are 
typically studied within Ericksen-Leslie theory, which as- 
sumes that the orientation of the normalized nematic di- 
rector is the only fluctuating variable This approxi- 
mation is adequate deep within the nematic phase, where 
the strength of nematic order is not significantly affected 
by thermal fluctuations. However, in the vicinity of 
the weakly flrst-order isotropic-nematic transition, sig- 
nificant fiuctuations in the nematic order are observed, 
suggesting that the phase-only approximation embodied 
in Leslie- Ericksen theory is inadequate . The study of 
nucleation in quenches from the isotropic to the nematic 
phase involves the growth of one phase within another, 
mandating the use of descriptions capable of describing 
both isotropic and nematic phases on the same foot- 
ing. In these and similar situations, a tensorial descrip- 
tion of nematic order which uses the symmetric, trace- 
less quadrupole moment tensor Qap, is appropriate, as 
first clarified by de Gennes in his Ginzburg-Landau the- 
ory of the isotropic-nematic transition'^. The Ginzburg- 
Landau-de Gennes (GLdG) approach provides a simple, 
but accurate phenomenological description of nematic 
fiuctuations in the static case[l|- 

The description of the fluctuating dynamics of the ori- 
entation tensor within GLdG theory has received consid- 
erably less attention 5, 6]. Understanding the results of 
inelastic scattering experiments on nematic systems [7[, 
the description of the rate of nucleation into the ne- 
matic phase the modeling of nontrivial stresses aris- 
ing from Casimir interactions Q and the calculation of 
the spectrum of capillary waves on the isotropic-nematic 
interface are all problems which require a dynami- 
cal theory of fluctuations in the orientation tensor. This 
problem is addressed in this paper, in which we present 
and solve the Langevin equations for dynamical fluctu- 



ations at equilibrium for the nematic orientation ten- 
sor. These are stochastic non-linear partial differential 
equations for the flve components of the orientation ten- 
sor. Analytical solutions can be obtained when these 
equations are linearized. For the solution of the gen- 
eral non-linear equations, we propose an efficient numer- 
ical method, based on a stochastic generalization of the 
method of lines. We compare our results with analytic 
results where such calculations are possible, finding ex- 
cellent agreement. 



II. FLUCTUATING NEMATODYNAMICS 

Orientational order in the nematic phase is described 
by a second-rank, symmetric traceless tensor Qafj{'x.,t). 
This is the second moment of the microscopic orienta- 
tional distribution function. The tensor can be expanded 
as 
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Qap = -zSinaUp - T,5ap) + -^Tllalp - maTnp). (1) 



The three principal axes of this tensor, obtained by di- 
agonalizing Qai3 in a local frame, specify the direction of 
nematic ordering n, the codirector 1 and the joint nor- 
mal to these, labeled by m. The principal values S and 
T represent the strength of ordering in the direction of n 
and m, quantifying, respectively, the degree of uniaxial 
and biaxial nematic order. 

The static fiuctuations of Qap can be calculated from 
a Ginzburg-Landau functional, first proposed by de 
Gennes, based on an expansion in rotationally invariant 
combinations of Qap and its gradients. The Ginzburg- 
Landau-de Gennes functional F is 

ic(TrQ2)2 (2) 
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Here, A = Ao{l - T/T*) T* denoting the supercoohng 
transition temperature, Aq a constant, Li is an elas- 
tic constant and a, I3,j denote the Cartesian directions. 
From the inequahty i(TrQ^)^ > (TrQ^)^, higher pow- 
ers of TrQ'^ can be excluded for the description of the 
uniaxial phase. Uniaxial phases are described by E' = 
while biaxial phases require E' ^ 0. For the nematic 
phase (rod-like molecules) B < Q whereas for the dis- 
cotic phase (plate-like molecules) B > 0. The quantities 
C and E' must always be positive to ensure bounded- 
ness and stability of the free energy in all phases. We 
omit other symmetry-allowed gradient terms in this pa- 
per, thus working in the limit where all three Frank con- 
stants are assumed to be equal. Such symmetry-allowed 
terms, as also total derivative surface terms, can be ac- 
counted for without essential change, using the numerical 
method described below. 

In the limit that hydrodynamic interactions may be 
neglected, i.e. the Rouse or free-draining limit, the dy- 
namical fluctuations of Qa/s are not coupled to other hy- 
drodynamic variables. The Langevin equations are those 
appropriate to a non-conserved order parameter with an 
overdamped, relaxational dynamics of the form 

5F 

dtQaP — —^aPfivjr: \-^al3, (3) 

Here the kinetic coefficients Fq.^^^, defined as Tapf^i, = 
^[Sa^iSpy + SaySp^ - \Sap5^u\-, cusurc that the dynam- 
ics preserves the symmetry and tracelessness property of 
the order parameter. In the absence of long-range forces, 
a local approximation for the kinetic coefficients is ad- 
equate and F can be taken as constant. The S^ap are 
symmetric, traceless Gaussian white noises, which sat- 
isfy a fluctuation-dissipation relation at equilibrium of 
the form 

(ec/3(x,t)) = 0, (4) 
(ea/3(x,t)e;..(x',i')> = 2kBTT^p^,J{^-^)5{t-t'\b) 

Here fc^ is the Boltzmann constant, T the temperature 
and denotes the average over the probability distri- 
bution of the noise. These Langevin equations, together 
with the fluctuation-dissipation relation for the noise, en- 
sure that the stationary one-point probability distribu- 
tion oiQap^ P[Qai3\, convergcs to Boltzmann equilibrium 
with P[Qc.p] - exp{-F/kBT). 

The equations above are five coupled, non-linear 
stochastic partial differential equations, with a noise term 
which has a tensorial structure. A numerical method 
of solution must maintain the symmetry and traceless 
of Qap. To ensure equilibrium dynamics, it must also 
maintain the balance between fluctuation and dissipa- 
tion. These two stringent requirements may be satisfled 
by transforming to a basis in which QajS is traceless and 
symmetric by construction. Symmetry and tracelessness 
of Qap is automatic. As we show below, the noise can be 
constructed out of independent Gaussian noises. 



We expand the orientational tensor in a basis of sym- 
metric traceless matrices T^^ as 

5 
i=l 

with Ti = y372'2¥,T2 = /T72(xx-y y),T3 = 

^2 'if , T'' = \/2 ¥¥ and = \/2 . The com- 
plete basis of matrices is orthogonal in the sense that 
T^^T^^ = Sij . In previous work we have presented explic- 
itly the equations for the basis coefficients (x, t) that 
follow from the deterministic part of the relaxational ki- 
netics dtQajB — ~^ai3tj.u5F/SQ^i, [ll| . (Thesc differ from 
the equations derived by others in that we include all 
non-linearities as well as an additional symmetry-allowed 
gradient (L2) term.) Here we focus on how an explicit 
construction of the noise can be implemented by expand- 
ing in the same basis. 
We expand the noise as 
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^„^(x,i)=^e.(x,07^:^, (7) 

i=l 

where each (x, t) is a zero- mean Gaussian white noise. 
From the orthogonality of the basis the inverse relation 
is 

^,(x,0=^?a/3(x,i)T^^. (8) 

a, 13 

From this, and the fiuctuation-dissipation relation it fol- 
lows that 

U^M,i^',t')) = 5](e„^(x,0C^.(x',t'))r:^7^^., (9) 

= 2kBTr6,jS{x- x')S{t - t'). 

This shows that the non-trivially correlated noise £,ai3 can 
be constructed from uncorrelated noises f^. Thus, by 
construction, the noise is symmetric, traceless and 
satisfies the fluctuation-dissipation relation. 

When anharmonic terms are ignored in the Ginzburg- 
Landau-de Gennes functional, the Langevin equations 
are linear and correlation functions can be calculated ex- 
plicitly in the T basis. Then, the Langevin equations of 
motion in terms of 

ai{q,t) = J d'^a;exp(— iq • x)a.i(x, i), (10) 

are 

dta,{q,t) ^ -T{A + L,q^)a,iq,t)+Uq,t)- (H) 

From this, the static and dynamic correlations follow im- 
mediately, 

C,j(q) = (a,(q)aj(-q)) = . ^ S,j, (12) 
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FIG. 1: (Color online) Autocorrelation function for the 
Ornstein-Uhlenbeck process, showing {v{t)v{t + T)) as a func- 
tion of the time increment r. The inset shows the histogram of 
fluctuations. It is Gaussian with the expected variance. The 
numerical parameters chosen are F = fcsT = 0.1, dt = 1.0. 
The average is taken over 10 independent realizations while 
the integration is performed for lO'^ SRK4 steps. 

and 

C»j(q,T) = (a,(q,i)aj(-q,i + r)) (13) 
= Cy(q)cxp[-r(yl + Liq2)r], 

where Cij(q) is the static structure factor. The static 
and dynamic correlations for Qap are then obtained by 
returning to the original basis. The stationary probabil- 
ity distribution generated by the Langevin dynamics is 
Gaussian with zero mean and variance Cij(q), consistent 
with Boltzmann equilibrium. 

III. STOCHASTIC METHOD OF LINES 

The fluctuating nematodynamics equations contained 
in Eq. ([3]) are five non-linear stochastic partial differen- 
tial equations. In general these have no analytical so- 
lutions and reliable numerical methods are therefore es- 
sential for their study. Here we combine the method of 
lines for solving initial- value partial differential equations 
with a stochastic Runge-Kutta integrator for systems of 
stochastic ordinary differential equations. This enables 
us to construct an accurate and efficient solver for the 
equations of fluctuating nematodynamics. Our results 
here build on previous work [Tlj . where a method of 
lines approach was used to solve the deterministic time- 
dependent Ginzburg-Landau equations numerically. The 
methodology here can thus be thought of as a general- 
ization of the method of lines to stochastic partial differ- 
ential equations. 

The method of lines is based on the idea of semidis- 
cretisation, where an initial-value partial differential 



equation in s pac e and time is discretised only in the spa- 
tial variable This yields a (possibly large) system 
of ordinary differential equations which is then solved by 
standard numerical integrators. To apply this method to 
stochastic partial differential equations, we must account 
for the fact that integrators for ordinary differential equa- 
tions do not automatically provide efficient and accurate 
solutions of stochastic differential equations. Qualita- 
tively, the noise term in a stochastic differential equation 
is a rapidly varying function and hence must be inte- 
grated with some care. At a more technical level, the 
noise is a Wiener process and the theory of stochastic 
integration must be used to evaluate it correctly (isj . 

Common stochastic integrators include those due to 
Maryuama [T^ and Milstein [is']. In this work, we use 
an integrator proposed recently by Wilkie [16], based on a 
multi-step Runge-Kutta strategy. The integrator is accu- 
rate and easy to implement by making small changes to 
a deterministic Runge-Kutta integrator. Further, since 
it is an explicit integrator, no matrix inversions are in- 
volved. This makes it attractive when the method of 
lines discretisation produces a large system of ordinary 
differential equations, as in our case. 

To test Wilkie's algorithm for a stochastic Runge- 
Kutta integrator, henceforth denoted as SRK4, we first 
check that the fluctuation-dissipation is obeyed. We 
performed a simple benchmark test on the Ornstein- 
Uhlenbeck process, represented as the Ito differential 
equation, 

dv{t) = -rv{t)dt + dW{t), (14) 

where dW{t) = y/Tt^TTdi Af (0,1) is the increment of 
the stochastic variable in the interval dt, and Af{0, 1) is 
a zero-mean unit- variance normal deviate. By construc- 
tion, the increments of this stochastic variable are inde- 
pendent and normally distributed with mean {dW{t)) = 
0. The particular choice of the variance ensures that the 
equilibrium distribution of is a Gaussian with variance 
fc^T. The stationary two-point autocorrelation of the 
velocity from Eqn. p^ is, 

{v{t)v{t + r)) = fcsTexp(-FT). (15) 

Fig.([T]) shows the autocorrelation as a function of time 
and the histogram of equal-time fluctuations of v. The 
variance (w^) = fc^T is correctly reproduced, as is the 
exponential decay of the autocorrelation function. We 
conclude that SRK4 is suitable as an integrator for prob- 
lems where the fluctuation-dissipation relation must be 
maintained. 



IV. NUMERICAL METHOD 

We now apply the method of lines together with 
SRK4 to obtain a stochastic method lines discretisa- 
tion (SMOL) for the equations of fluctuating nemato- 
dynamics. We benchmark our numerical results by com- 
paring autocorrelations within a harmonic theory which 
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FIG. 2: (Color online) Histogram of the real part of ai(q) for 
Tlx = ny = 6 in a. box of dimension Lx = Ly — 16, using a 
harmonic free energy, with ksT = A — 0.05 and Li — 0.5, 
r = 1.0. The fluctuations are Gaussian with the expected 
variance. The histogram is obtained from 20 independent 
realizations, each realization contributing 4000 time steps. 



accurately describes fluctuations about the isotropic 
phase. We then consider expansions about the ordered 
state, comparing static correlations obtained analytically 
within the Frank approximation with our numerical re- 
sults. 

We use a finite-difference discretisation with nearest- 
neighbour stencils for gradients and the Laplacian. We 
implement periodic boundary conditions. Specifically, in 
three dimensions, we consider a box of dimension L^, Ly 
and Lz along the Cartesian directions, and grid these 
lengths with equal grid spacing Aa; = Ay — Az — 1. 
The latter defines lattice units for the spatial coordi- 
nate. We define corresponding discrete time units for 
the temporal variables by choosing At = 1. Fourier 
modes are labelled by the wave- vector q = (fe, Qz), 
where each component is of the form qa — 2Tma/La, 
with ria = 0, 1, 2, ... , (Lq, — 1). 

With this discretisation, the Laplacian in Fourier space 
is given by 



£(q) = 2[cos(gx) + cos{qy) + cos{qz) - 3]. 



(16) 



The nearest-neighbour finite difference stencil suffers 
from lack of isotropy at high wavenumbers. This can be 
improved through the use of higher-point stencils [l"7l - [l9| . 

Applying the method of lines discretisation to Eq. ^ 
reduces it to a system of stochastic ordinary differential 
equations, whose Fourier representation in the harmonic 
approximation of Eq. (jlip is 



(17) 



The Fourier representation of the drift-diffusion dynam- 
ics is encoded in the linear operator V, 
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FIG. 3: (Color online) (a) Contour plot of the structure factor 
Cij(q) and (b) angular average of Ci-, (q). The parameters are 
ksT = A ^ 0.05, Li = 0.5, r = 1.0. The time averaging is 
over 5 x 10* time steps and ensemble averaging is over 20 
independent realizations in a box with = Ly = 64. The 
relaxation time scale is r = {rA)~^ — 20, the diffusion time 
scale of the smallest Fourier mode is Td = Z/^/(47r^Lir) = 
207.51 and the correlation length is A = ^Li/A = 3.16. 



Fourier representations of the one and two-dimensional 
method of lines discretisations are obtained by setting 
the corresponding wavenumbers to zero. The static and 
dynamic autocorrelations in Fourier space follow in a 
straightforward manner though the replacement of q^ by 
its discrete Laplacian representation. The results are 



a,(q) 



Qj(q)exp(-r2?T). 



(19) 
(20) 



V = A-LxC{q). 



(18) 



It is also useful to define an angle-averaged structure fac- 
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q) for comparison with the nu- 



FIG. 4: (Color online) Autocorrelation Cij(q, r) for the linear 
Langevin equation, calculated for small wavenumbers. Nu- 
merical parameters are fcsT = A — 0.05, Li = 0.5, F = 1.0. 
The integration is performed for 4.2 x 10'^ time steps on a 16^ 
grid. The time average is taken over 4 x 10^ time steps and 
an ensemble average is taken over 40 independent realizations. 
The relaxation time scale is r = {rA)~^ = 20, the diffusion 
time scale of the shortest Fourier mode Td — L^/(47r^LiF) ^ 
13 and the correlation length A — ^ L\ /A — 3.16. 
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FIG. 5: (Color online) Angularly averaged structure factor 
of the nematic phase. The numerical parameters are fcsT = 
0.05, A = -3.5, B = -10fcsr,C = 2.67, Li = 32.0, F = 0.01 
on a 64^ grid. The relaxation time scale is r = {rA)~^ = 
28.57, the diffusion time scale of the shortest Fourier mode is 
Td = L^/(47r^LiF) — 324.23 and the correlation length is A = 
■yjLr/A = 3.02. The time average is taken over 5 x 10* time 
steps and an ensemble average is taken over 20 independent 
realizations. 



merical simulation. 

We now compare theoretical and numerical results: In 
Fig. ^ we show the histogram of the ai{q) for a par- 
ticular Fourier mode. This is normally distributed, as 
expected, with zero mean and variance as required by 
thermal equilibrium. Similarly, all Fourier modes exam- 
ined have correct normal distributions. The variances 



obtained are compared in Fig. 3(a) with the analytical 
values by plotting contours of Cij(q). There is excel- 
lent agreement. A close inspection reveals some degree 
of anisotropy in both the analytical and numerical results 
at high wavenumbers. This is attributed to the lack of 
isotropy of the nearest-neighbour finite-difference Lapla- 
cian mentioned earlier. However, the anisotropics are 
removed upon angular averaging, as shown in Fig. |3(b)[ 
Thus, the present discretisation should be adequate in 
most cases, unless highly accurate isotropics are required 
from the simulation. From these results, we conclude 
that correlations in thermal equilibrium are accurately 
captured by the stochastic method of lines approach. 

We next compare the dynamics of fluctuations at equi- 
librium, by comparing two-point autocorrelation func- 
tions calculated analytically and numerically. Fig.(|4]) 
shows Cy (q, r) for three sets of Fourier modes. The ex- 
ponential decay of the autocorrelation function is repro- 
duced accurately within the numerics and fit the theo- 
retical curve very closely. We conclude, therefore, that 
the stochastic method of lines accurately reproduces both 
static and dynamic fluctuations in a harmonic theory. 

Finally, we compare theory and simulation in a situ- 
ation where a linearization of the Qa/3 equations about 
Qaf3 = is inapplicable, that of director fluctuations 
within the nematic phase. In the T basis, the equations 
of motion are 



dta, = -T[{A + CTrQ^)a, 



(21) 



where, is the traceless symmetric projection of 

Q^^. These equations of motion are anharmonic, and the 
analytical solutions obtained earlier within the harmonic 
expansion are no longer available for comparison. We 
therefore extract the fluctuations of the angular displace- 
ments from the uniform nematic ground state using an 
approach based on the Frank free energy. 

Consider an uniform uniaxial nematic with director no 
with small fluctuations (5n(x). Decomposing the fluctu- 
ations into parts parallel and perpendicular to no and 
imposing the normalization of the director, we find from 
the Frank free energy that, 



(l'5n^(q)P) = 



(22) 



where K = {9S^ /2)Li is a Frank constant [4]. Since fluc- 
tuations in the plane perpendicular to ng can be charac- 
terized through a single angle 9, an equivalent result is 
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(6l(q)6'(-q)) = ksT/Kq^. In the semidiscrete represen- 
tation, we obtain 

Fig-® shows the angular average of the static 
correlations of director fluctuations Gg{q) — 
^lql^^(6'(q)0(— q)). The formally divergent q = 
mode is excluded both from the numerical data and 
analytical result. Given that the analytical result is 
obtained from a linearization about the aligned state 
whereas the numerical solution is calculated from the 
equations of motion arising from the full non-linear free 
energy, the agreement between theory and simulation 
is satisfactory. The stochastic method of lines thus 
accurately captures equilibrium fluctuations in the 
ordered state as well as in the disordered one. 



V. DISCUSSION AND CONCLUSION 

The numerical method of solution presented in this 
paper can be applied to a variety of problems in ne- 
matodynamics where accounting for fluctuations in ne- 
matic order are important. To the best of our knowledge, 
no systematic study of anharmonic fluctuations exists 



within the time-dependent Ginzburg-Landau framework. 
In previous work, Stratonovich Q presented fluctuating 
equations of motion for harmonic fluctuations in terms 
of Langevin equations, analyzing these equations within 
the traceless symmetric basis described here in a straight- 
forward manner Q. Our equations of motion contain the 
necessary nonlinearities and our numerical methodology 
accounts for them in a computationally straightforward 
way. 

We also present, to the best of our knowledge for the 
first time, a systematic stochastic integration scheme ca- 
pable of yielding highly accurate solutions of the non- 
linear equations of nematodynamics. These equations 
can be applied to study fluctuations of nematic order at 
the isotropic-nematic interface, pseudo-Casimir interac- 
tions close to the isotropic-nematic transition, as well as 
nuclcation phenomena in nematogens. The study of these 
and similar problems is ongoing and will be reported else- 
where. 
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